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Abstract 

The two-dimensional harmonic polylogarithms G{a{z);y), a generalization of the harmonic poly- 
logarithms, themselves a generalization of Nielsen's polylogarithms, appear in analytic calculations of 
multi-loop radiative corrections in quantum field theory. We present an algorithm for the numerical 
evaluation of two-dimensional harmonic polylogarithms, with the two arguments y, z varying in the 
triangle 0<?/<l, O<0<1, Q < {y + z) < 1. This algorithm is implemented into a FORTRAN 
subroutine tdhpl to compute two-dimensional harmonic polylogarithms up to weight 4. 



PROGRAM SUMMARY 



Title of program: tdhpl 
Version: 1.0 
Release: 1 
Catalogue number: 

Program obtained from: Thomas. GehrmannOcern.ch, Ettore.RemiddiQbo.infn.it 
E-m,ail: Thomas . GehrmannOcern. ch, Ettore.RemiddiObo.infn.it 

Licensing provisions: no 

Computers: all 

Operating system: all 

Program language: F0RTRAN77 

Memory required to execute: Size; 1144k 

No. of bits per word: up to 32 

No. of lines in distributed program: 34589 

Other programs called: hplog (from the same authors, Comput. Phys. Commun. 141 (2001) 296) 

External files needed: none 

Keywords: harmonic polylogarithms, Feynman integrals 

Nature of the physical problem: numerical evaluation of the two-dimensional harmonic polylogarithms up 
to weight 4 for real arguments restricted to < y < 1, < ^; < (1 — y). These functions are emerging in 

Feynman graph integrals involving three mass scales. 

Method of solution: for small values of the argument: series representation with expansion coefficients 
depending on the second argument; other values of the argument: transformation formulae. 

Restrictions on complexity of the problem: limited to 2dHPLs of up to weight 4, the algorithms used here 
can be extended to higher weights without further modification. 

Typical running time: on average 1.8 ms (both arguments varied)/0.5 ms (only one argument varied) for 
the evaluation of all two-dimensional harmonic polylogarithms up to weight 4 on a Pentium III/600 MHz 
Linux PC. 
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LONG WRITE-UP 



1 Introduction 

The two-dimensional harmonic polylogarithms, a generaUzation of the harmonic polylogarithms H(a, x), 
have been introduced in for the analytic evaluation of a class of two-loop, off-mass-shell scattering 
Feynman graphs in massless QCD. 

Two-dimensional harmonic polylogarithms are obtained by the repeated indefinite integration of rational 
fractions, involving one further independent variable besides the integration variable. The first appearance 
of functions of this type in the mathematical literature was in a series of articles of E.E. Kummer |^ in 
1840. Allowing for arbitrarily many different variables to appear in the rational fractions, one obtains a 
class of functions introduced by Poincare, which are called 'hyperlogarithms'. These have been studied in 
great detail in the works of J. A. Lappo-Danilevsky ||^. Hyperlogarithms and their generalization 'multiple 
polylogarithms' are receiving renewed attention in the mathematical literature (see for instance the review 
of A.B. Goncharov [Q). 

It has been known for a long time that the analytic evaluation of integrals in perturbative quantum 
field theory gives rise to the Euler dilogarithm Li2(a;) and its generalizations, Nielsen's polylogarithms [^. 
A reliable and widely used numerical representation of these functions (GPLOG) |^ has been available for 
already thirty years. Going to higher orders in perturbation theory, it was recently realized that Nielsen's 
polylogarithms are insufficient to evaluate all integrals appearing in Feynman graphs at two loops and 
beyond. This limitation can only be overcome by the introduction of further generalizations of Nielsen's 
polylogarithms. This generalization is made by the harmonic polylogarithms (HPLs) |Q, appearing in loop 
integrals involving two mass scales, and two-dimensional harmonic polylogarithms (2dHPLs) [Q, appearing 
in loop integrals involving three mass scales. These functions are already now playing a central role 
in the analytic evaluation of Feynman graph integrals [^|^]. HPLs appear moreover as inverse Mellin 
transformations of harmonic sums, which were investigated and implemented numerically in [^, while 
2dHPLs and multiple polylogarithms also appear for example if generalized hypergeometric functions are 
expanded in their indices around integer values |p^ . 

A subroutine (hplog) for the numerical evaluation of HPLs for arbitrary, real values of the argument 
was presented in 

In this paper, as a continuation of we briefly review the analytical properties of the 2dHPLs and 
then show how those properties can be used for writing a FORTRAN code that evaluates the 2dHPLs up to 
weight 4 (see Section below for the definition of the weight; 4 is the maximum weight required in the 
calculations of [Q) with absolute precision better than 3 x 10~^^ {i.e. standard double precision) with a 
few dozens of elementary arithmetic operations per function. Given the large number (256) of 2dHPLs of 
weight 4, the many algebraic relations among them, and the fact that any application is likely to involve a 
large number of them at the same time (see for instance the results of [Q) our FORTRAN routine evaluates 
the whole set of 2dHPLs up to the required weight - as in hplog [|ll|, but at variance with GPLOG ||], which 
evaluates separately (and up to weight 5) the various Nielsen polylogarithms. While hplog evaluated the 
HPLs for arbitrary real value of the argument, yielding complex results, we restrict the arguments of the 
2dHPLs to the region 0<2/<l — 2:,0<z<l, where the 2dHPLs are real. This region is the only relevant 
region for applications in quantum field theory ||l| . 

The plan of the paper is as follows. Section || recalls the definitions of the HPLs. Their algebraic 
properties are discussed in Section |^, where we show how to use these properties for separating the functions 
into reducible and irreducible ones. In Section ^, we discuss the behaviour of 2dHPLs for special values of 
the argument. Relations between 2dHPLs for different ranges of the arguments are derived in Section |^. 
Section ^ studies the analytic properties that allow the performing of converging power series expansions 
and the acceleration of their convergence. Section ^ explains how the properties recalled above are used to 
implement the 2dHPLs into the FORTRAN subroutine tdhpl, and Section ||how the correct implementation 
is checked. Finally, we describe the usage of the subroutine tdhpl in Section and provide a few numerical 
examples in Section |l^. We enclose two appendices. Appendix ^ compares the notations used for 2dHPLs 
in different previous publications and Appendix |^ provides relations between particular cases of the 2dHPLs 
and Nielsen's generalized polylogarithms. 
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2 Definitions 



The 2dHPLs family which we consider here is obtained by the repeated integration, in the variable y, of 
rational factors chosen, in any order, from the set 1/y, l/{y — l), l/(y + z — 1), l/(y + z), where z is another 
independent variable (hence the 'two-dimensional' in the name). It is clear that the set of rational factors 
might be further extended or modified; for the harmonic polylogarithms discussed in |pl]| the set of rational 
factors was for instance 1/y, !/(?/ — 1), 1/iy + 1), involving only constants and no other variable besides y. 
More precisely and in full generality, let us define the rational factor as 

g(a;y) = ^, (2.1) 
y-a 

where a is the index, which can depend on z, a = a{z); the rational factors which we consider for the 
2dHPLs are then 



g(0;y) 
g(i;y) 

;(1 - z;y) 



1 

y 



1 



y-1 
1 



y + z - 1 



g(-;y) - (2.2) 

With the above definitions the index takes one of the values 0, 1, (1 — z) and {—z). 

Correspondingly, the 2dHPLs at weight w — 1 (i.e. depending, besides the variable y, on a single further 
argument, or index) are defined to be 

G{0;y) = Iny, 
G(l;y) = ln(l-j/), 

G{1 -z;y) = ln(^l-^ 

G{-z;y) = ln(l + ^) . (2.3) 

The 2dHPLs of weight w larger than 1 depend on a set of w indices, which can be grouped into a w- 
dimensional vector of indices a. By writing the vector as a = (a, b), where a is the leftmost component of a 
and b stands for the vector of the remaining {w ~ 1) components, the 2dHPLs are then defined as follows: 
if all the w components of a take the value 0, a is written as 0^ and 

G(0^;2/) = ^ln'"y, (2.4) 



y 



G{d;y)^ / dy'g{a;y') G{b;y') . (2.5) 



while, if a 7^ Oui, 



In any case the derivatives can be written in the compact form 

^G{a;y) = gia;y)Gib;y) , (2.6) 

where, again, a is the leftmost component of a and b stands for the remaining (w — 1) components. 

From ( ^.4| ) and (|]^) , one arrives immediately at a multiple (or repeated) integral representation of the 
2dHPL: 

G(m„;y)=/ dtig{mi;ti) dt2g{m2;t2) ■ . ■ dt^g{m^;t^) , (2.7) 
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valid for rriu, 7^ 0, and 

G(m.u,;y)= / dtig{mi;ti) dt2g(m2; ^2) ■ • • / dtyg{mv]tv)G{Ow-v;tv) , (2.8) 
Jo Jo Jo 

valid for ifiw = [iTiv, Ow-v) with rriv 7^ O^,. 



The definition is essentially the same as for the harmonic polylogarithms of 1 1 1 1 , if allowance is made 
for the greater generality of the 'indices', which can now depend on the second variable z. Let us, however, 
stress an important difference between the present definitions and the notation already used in [Q , where 
the rational factors were indicated by i{a,x) and the harmonic polylogarithms by H(a, a;); we have indeed 

f(l;x) = -g(l;a;), 
f(l - z;x) = -g(l - z;x) , 
= g{-z;x) , 

(2.9) 



while there is no change when a — 0: 



i{0;x) = g(0;x). (2.10) 



Also for a = —1 one would have f(— — g(— l;a;), but we will not consider this case here as it never 
appears together with the other values of the indices (1 — z), {—z) in The same applies between the 
harmonic polylogarithms H previously introduced ^ and the 2dHPLs, as any H-function of |l] goes into 
the corresponding G-function, with the following correspondence rules: the indices (1), (1 — z) of H remain 
unchanged as indices of G, but each occurrence of (1), (1 — z) gives a change of sign between H and G; any 
index (z) of H goes into an index (— z) of G (which keeps the same sign as H). One has for instance 

H(z,l-z;y) = -G(-z, 1 - z; y) , 
H(0,z,l-z,l;y) = G(0, -z, 1 - z, 1; y) , (2.11) 

and so on. The main advantage of the new notation is the (obvious) continuity in z of the g's and the G's; 
one has for instance 

limg(l-z;y) = g(0;y) , (2.12) 

to be compared with 

limf(l-z;y) = -f(0;y) , (2.13) 

and the same applies to any index of a G-function (when the limit exists). Note, however, that the positivity 
for positive value of the argument is lost - so that, for instance, one has G(0, 1; 1) — — 7r^/6, to be compared 
with the more elegant relation H(0, 1; 1) = 7r^/6. 

The 2dHPLs can also be viewed as a special case of the hyperlogarithms, which are discussed frequently 
in the mathematical literature. We summarize the various available notational conventions in Appendix A, 
where we also provide appropriate translation rules. 



3 The algebra and the reduction equations 

Algebra and reduction equations of the 2dHPLs are the same as for the ordinary HPLs. They were discussed 
extensively in 0,0, and identical formulae apply regardless of the actual values of the indices. In the 
following, we briefly summarize the results of |0,|ll[l, without providing explicit examples. 

The product of two 2dHPLs of a same argument x and weights p, q can be expressed as a combination 
of 2dHPLs of that argument and weight r — p + q, according to the product identity 

G{p;x)G{q;x) = ^ G(f;a;), (3.1) 
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where p, q stand for the p and q components of the indices of the two 2dHPLs, while p l±) g represents all 
possible mergers of p and (f into the vector r with r components, in which the relative orders of the elements 
of p and q are preserved. 



Another class of identities is obtained by integrating (2.4) by parts. These integration-by-parts (IBP) 
identities read: 

G(mi, . . . , TO^; a;) — G(toi; a;)G(m2, . . . , m^; x) — G(m2, mi; a;)G(m3, . . . , m^; x) 

+ ... + (-l)«+iG(TO„...,mi;x) . (3.2) 

These identities are not fully linearly independent from the product identities. 



By using Eq. (3J) at weight w = 2 for all possible independent values of the indices, one obtains 
10 independent relations between the 16 2dHPLs of weight 2 and products of 2 2dHPLs of weight 1; 
those relations can be used for expressing 10 of the 2dHPLs of weight 2 in terms of 6 2dHPLs of weight 
2 and products of 2 2dHPLs of weight 1. The choice of the 6 2dHPLs (referred to, in this context, 
as 'irreducible') is by no means unique; by choosing as irreducible 2dHPLs of weight 2 the 6 functions 
G(0, l;j/), G(0, 1 — z;2/), G{0,—z;y), G(l — z,l;y), G{—z,l;y), G(—z, 1 — z; y), the reduction equations 
expressing the 10 'reducible' 2dHPLs of weight 2 in terms of the irreducible 2dHPLs read 

G(-z,0;2/) = G(0;y)G(-z;y)-G(0,-z;2/), 
G(l-z,0;y) = G(0; y)G(l - z; y) - G(0, 1 - z; y) , 
G(l,0;y) = G(0;y)G(l;y)-G(0,l;y), 

G(-z,-z;y) = ^Gi-z;y)G{-z;y) , 

G{1 -z,~z;y) = G(l - z; y)G(-z; y) - G(-z, 1 - z; y) , 
G(l,-z;2/) = G(l;y)G(-z;y)-G(-z,l;y), 

G(l-z,l-z;y) = ^G{1 - z;y)G{l - z:y) , 

G(l,l-z;y) = G(l;y)G(l-z;y)-G(l-z,l;y), 

G(l,l;y) = iG(l;y)G(l;y), 

G(0,0;y) = ^G{0;y)G{0;y) . (3.3) 

Similarly, at weight 3, one has 64 2dHPLs and 44 independent product and IBP identities, expressing 44 
reducible 2dHPLs in terms of 20 irreducible ones; at weight 4 there are 256 2dIIPLs, 196 independent 
identities, and correspondingly 196 reducible and 60 irreducible 2dHPLs. 

The 2dHPLs chosen as the irreducible set in the program described here are listed in Table |l|; for 
convenience of later use (see in particular Section ^ we grouped them according to the different combina- 
tions of the occurring indices. Note in particular that, with our choice, the index (0) is never present as 
rightmost (or trailing) index of the 2dHPLs of the irreducible set (except for the trivial case, at w = 1, of 
H(0;y) = lny). 



4 Special values of the argument 

Aty = l — z, y = —z and y = 1, the 2dHPLs can be expressed in terms of HPLs of argument z. 
To obtain such expressions, one can write the obvious equation 

G(m(z);y(z))-G(rn(z = 0);y(z = 0))-f / dz'/;G(m(z'); y(z')) , (4.1) 

Jo 

where y(z) stands for any of the above particular values y = 1 — z, y = — z or y = 1 (the argument applies 
for y taking the constant value 1 as well), m(z) is any set of indices and it is understood that the boundary 
z = is to be replaced by z = 1 if G(m(z — 0);y(z = 0)) is divergent. The derivative d/dz' is then 
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Weight 


Indices 




w = 1 




G(0; y); G(l; y); G(l - 2; y); G(-2; y) 


w = 2 


(0,1) 


G(0, l;y) 




(0, 1 - z) 


G(0, 1 - 2; J/) 




(0, -z) 


G(0, -2; J/) 




(0, 1 - z, 1) 


G(l - 2, l;y) 




(0, -z, 1) 


G(-2, l;i/) 




(0, —z, 1 — z) 


G(-2, 1 - 2;y) 


w = 3 


/n 1 \ 

(0- 1) 


G(0, 0, l;j/); G(0, 1, l;j/) 




(0, 1 - z) 


G(0, 0, 1 — 2; y); G(0, 1 — 2, 1 — 2; J/) 




(0, -z) 


G(0, 0, — 2; J/); G(0, — 2, — 2; J/) 




(0, 1 — z, 1) 


G(0,1 — 2,1;?/); G(l — 2,0,l;y); G(l — 2,1,1; y); G(l — 2,1 — 2,1;?/) 




(0, -2, 1) 


G(U, — 2, 1; y); o(— 2, 0, 1; y); G(— 2, 1, 1; ?/); 0(— 2, — 2, 1; y) 




(0, -z,l - z) 


G(0, -2, 1 - 2;y); G(-2, 0, 1 - 2;y); G(-2, 1 - 2, 1 - 2;y); G(-2, -2, 1 - 2;y) 




(0, -2, 1 - 2, 1) 


G(l-2,-2,l;y); G(-2, 1 - 2, 1; y) 


w = A 


(0,1) 


G(0,0,0,l;y); G(0, 0, 1, 1; y); G(0,l,l,l;y) 




(0, 1 - 2) 


G(0, 0, 0, 1 - 2; y); G(0, 0, 1 - 2, 1 - 2; y); G(0, 1 - 2, 1 - 2, 1 - 2; y) 




(0, -2) 


G(0, 0, 0, -2; y); G(0, 0, -2, -2; y); G(0, -2, -2, -2; y) 




(0,1 -2,1) 


G(0, 0, 1 - 2,l;y); G(0, 1,1-2, l;y); G(0, 1-2,0, l;y); 
G(0, 1-2,1, l;y); G(0, 1 - 2, 1 - 2,l;y); G(l - 2,0,0, l;y); 
G(l - 2,0, 1, l;y); G(l - z,0, 1 - 2,l;y); G(l - 2,1,1, l;y); 
G(l - 2,1 - 2,0, l;y); G(l - 2, 1 - 2,1, l;y); G(l - 2, 1 - 2, 1 - 2, l;y) 




(0,-2,1) 


G(0, 0, -2, 1; y); G(0, 1, -2, 1; y); G(0, -2, 0, 1; y); 
G(0, -2, 1, l;y); G(0, -2, -2, 1; y); G(-2, 0, 0, 1; y); 
G(-2,0,l,l;y); G(-2, 0, -2, 1; y); G(-2, 1, 1, 1; y); 
G(-2, -2,0, l;y); G(-2, -2, 1, 1; y); G(-2, -2, -2, 1; y) 




(0,-2,1-2) 


G(0,0, -2, 1 - 2;y); G(0, 1 - 2, -2, 1 - 2;y); G(0, -2,0, 1 - 2;y); 
G(0, -2, 1 - 2,1 - 2;y); G(0, -2, -2, 1 - 2;y); G(-2, 0, 0, 1 - z;y); 
G(-2,0, 1 - 2,1 - 2;y); G(-2, 0, -2, 1 - z;y); G(-2, 1 - 2, 1 - 2, 1 - 2;y); 
G( — 2,— 2,0, 1 — 2;y); G( — 2,— 2,1 — 2,1 — 2;y); G( — 2, —2, —2,1 — 2; y) 




(0,-2,1-2,1) 


G(0, 1-2, -2, 1; y); G(0, -2, 1 - 2, 1; y); G(l - 2, 0, -2, 1; y); 
G(l - 2, 1 - 2, -2, l;y); G(l - 2, -2,0, l;y); G(l - 2, -2, 1, l;y); 
G(l - 2, -2, 1 - 2, 1; y); G(l - 2, -2, -2, 1; y); G(-2, 0, 1 - 2, 1; y); 

G(-2, 0, -2, 1; y); G(-2, 1,1-2, 1; y); G(-2, 1 - 2, 0, 1; y); 
G(-2,l - 2, 1, l;y); G(-2, 1 - 2, 1 - 2,l;y); G(-2, 1 - 2, -2, l;y); 
G(-2, -2,1-2, l;y) 



Table 1: List of irreducible 2dHPLs chosen in the numerical implementation. 



carried out on the representation of G{m{z'); yiz')) as a repeated integral ( p.7| ) or ( [2.81 ). Quadratic factors 
of the form [g{mi;tj)]'^ , which can appear when carrying out the 2;'-derivative explicitly, are reduced to 
single powers by integration by parts and partial fractioning (to be iterated when needed). The resulting 
repeated integrals can then be identified as a combination of multiple integral representations of HPLs of 
argument z. 

As an example, we evaluate G(l, 1 — z; y) in y — 1 — z: 

G(l,l-z;l-z) = G(1,0;0)+ / dz'— G(l, 1 - z'; 1 - z') 

Ji dz 
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= H(l,0;z)-H(l,0;l) + H(0,0;z) 



(4.2) 



Note the appearance, in the above calculation, of H(l; 1 — z'), which is then replaced by — -^(0; z') (according 
to the very definition, see for instance (A.l)). In more complicated cases the replacement is carried out in a 
recursive manner, i.e. by using in the transformation of 2dHPLs of weight w the results obtained previously 
for 2dHPLs of weight less than w. 

The above sketched algorithm has been programmed in FORM [ p2[ to derive the transformation formulae 
in y = 1 — z, y — —z and y = 1 for all the irreducible 2dHPLs up to weight 4. 



5 Identities for related arguments 

In the context of the numerical evaluation of the 2dIIPLs, it is sometimes convenient to map G(to(z); y) into 
G{rn(z); y') with y' = 1 — z — y. Indeed, this mapping allows us to rewrite a 2dHPL whose argument y lies 
in the range (1 — z)/2 < y < [1 — z) in terms of 2dIIPLs with argument y' in the range < y' < (1 — z)/2. 
As a consequence, power series expansions for the 2dHPLs are needed only in the region 0<y<(l — z)/2, 
thus avoiding potentially non-analytic points at y = (1 — z). 

Much as in the previous section, the mapping can be obtained by writing the obvious equation 

rv 



G{m{z); 1 - y - z) = G(m(z); 1 - z) + 



dy'AG(m(z);l-y' 
dy' 



(5.1) 



where the boundary y = can be replaced by y = 1 — z if G(m(z); 1 — z) is divergent (as is the case 
when 1 — z appears as leftmost index in m(z)). The y'-derivative is again carried out on the multiple 
integral representation of G(m(z); 1 — y — z) (2. 7), (2. 8). A repeated apphcation of integration by parts and 
partial fractioning then generates an expression that can be identified as a bilinear combination of 2dIIPLs 
G(m(z);y) and ordinary HPL H(m;z). As an example, consider 



G(l,l-z;l-y-z) 



G(l,l-z;l 
G(l,l-z;l 

G(l,l -z;l 



z) + 
z) + 



dy'— G(l,l-z;l 



y 



y d 
dy' 



y 



dy'- 



-1 



G(l,l-z;l-z) + 



-y - z 

y dy' 



h - 



1 



dt2 



t2 



dt2 







t2 + Z - 1 

[G(0;y')+H(l;z)] 



/o y' + z 

H(l, 0; z) - H(0, 1; 1) + H(0, 0; z) + G(-z, 0; y) + G(-z; y)H(l; z) 



(5.2) 



Note the use, in the last step, of the results of the previous section. Interchange formulae y ^ 1 — y — z have 
been derived for 2dHPLs up to weight 4 by programming the above sketched algorithm in FORM 12 1. As 
the transformation algorithm presented above, the interchange algorithm also works recursively, by using 
the interchange formulae obtained previously for lower weights. 



6 The analyticity properties 



Let us recall that we are interested only in the region < y < 1, 0<z<l, < y + z < 1. There are 
two possible right cuts in y, 

• y ~ 1, coming from g(l; y) — l/{y — 1) (the same as for the HPLs), 

• y = I — z, coming from g(l — z;y) — l/{y + z — 1), 
and a single left cut 

• y ^ —2, coming from g{—z; y) = + z) (remember z > 0) 

We use a set of irreducible 2dHPLs in which the indices appear in the order (0, {—z), (1 — z), 1), so that 
a rightmost (—z) index can have at its left only (0)'s or {—z)^s, a rightmost (1 — z) can have (0)'s, (— ^)'s or 
(1 — 2:)'s, and finally a rightmost (1) can have at its left (0)'s, {—z)^s, (1 — z)'s or (l)'s. As a consequence, 
a rightmost (0) and the related logarithmic cut at y = are never present in the 2dHPLs of the irreducible 
set chosen here (except in the trivial case, at w = 1, of H(0;y) = Iny). 

Following g, each 2dHPL can be written as 

Gim{z)-y) = G+(77i(z);y) + G_(m(z);y) , (6.1) 

where G+{'m{z);y) contains only one or both right cuts, and G_(m(z);y) contains only the left cut. Note 
that we keep together the two right cuts, when both are present. 

The separation of the cuts is carried out iteratively. At weight w = 1 the cut-structure is 

G(l;y) = G+(l;y), 
G(l -z;y) = G+(l-z;y), 

G(-z;y) - G_(-z;y), (6.2) 

as G_(l;y) = G_(l - z;y) = G+(-z;y) = 0. 

At higher weights, the following separation formulae apply: 



with 



G{a{z),rn{z);y) = G+{a{z),rn{z):y) + G-{a{z),m{z):y) , (6.3) 



r^G±(m(z);y'), 
Jo V 

y Ay' 



G± 


(0,m(z);y) = 


G+ 


(l,m(z);y) = 


G_ 


(l,m(z);y) = 


G+(l- 


-z,rh{z);y) 


G_(l - 


^z,rh{z);y) = 


G+(- 


-z,Tn{z);y) = 


G-(- 


-z,m{z);y) = 



y'-i 

y dy' 



(G+(m(z);y') + G_(m(z);l)) , 
I -;i^(G_(m(z);y')-G_(77i(z);l)) , 

"^y' (G+(m(z);y') + G_(m(z);l-z)) 
(G_(77i(z);y')-G_(m(z);l-z)) 



y'-i + z 
' dy- 

y'-i + z 
" Ay' 
y' + ^ 

" dy- 

y' + z 



(G+(m(z);y')-G+(77i(z);-z)) , 

(G_(77i(z);y') + G+(m(z);-z)) . (6.4) 



The cut corresponding to the rightmost index, see (5^), is called the primary cut, the other cut (when 
present) the secondary cut. The separation of cuts docs require the computation of the 2dHPLs at the 
points y = 1 — z, y = —z or y = 1, which was explained in Section 0. 
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It is important to note that, for a given 2dHPL G(to(z); y), the secondary cut consists of a product of 
H(a; z) times only 2dHPLs of weight fower than G{m(z);y). If the 2dIIPLs are evaluated for increasing 
weight, the secondary cut contribution does not need to be evaluated anew. Only the principal cut contains 
2dHPLs of the same weight as G(m(z); y), and has to be evaluated. 



As an example of the separation of the cuts according to the rules of (6^), let us consider 

G(0, -z, 1, 1; y) = G+(0, -z, 1, 1; y) + G_(0, -z, 1, 1; y) , (6.5) 

with 

G+(0,-z,l,l;y) - G(0, -z, 1, 1, y) - H(-l, -1; z)G(0, -z; y) , 

G„(0,-z,l,l;y) = H(-l, -1, z)G(0, -z, y) . (6.6) 



7 The numerical evaluation 

As in we evaluate the irreducible 2dHPLs first, and then the reducible ones by using the formulae 
expressing them in terms of the irreducible ones, see Section ||. 

For the purpose of the numerical evaluation of the 2dIIPLs, we restrict ourselves to the following values 
of the arguments, depicted by the outer triangle in Fig. |l|: 

0<y<l, 0<z<l, with y + z < 1 , (7.1) 

which is the only kinematical region relevant to the calculation of physical amplitudes in quantum field 
theory (note that there are other kinematically allowed regions in the (y, z)-plane, which can however be 
related to the above triangle by analytic continuation, see Appendix of 

When the vector of the indices of the 2dHPL contains only a single index different from (0), which can 
however be repeated several times, one can rescale the argument, obtaining in that way an HPL of suitable 
argument, which can be evaluated with the routine hplog [|lT| . 

If the argument y is in the inner triangle on the right of Fig. |l|, (1 — z)/2 < y < (1 — z), we use 
the results of Section || to express the 2dHPLs in terms of 2dHPLs of argument y' = (1 — y — z), with 
< y' < (1 -z)/2. 

Finally, when the argument y is in the triangle < y < (1 — z)/2 , on the left of Fig. |l|, we evaluate the 
primary cuts of the 2dHPLs of the irreducible set by their series expansion in powers of y around y = 0, 



their secondary cuts by their expression in terms of 2dHPLs of lower weight, see (6.3) and (|6^), and then 
sum the two cuts. 




Figure 1: Kinematic regions for the evaluation of the 2dHPLs 
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The coefficients of tlie expansions in y of the primary cuts are in turn functions of z. For Q < z < 1/2 
it is convenient to expand the coefhcients, which often show nominal divergences in 1/z, generating serious 
numerical instabilities, in powers of z, obtaining in that way a stable and quickly convergent expression. 
For some of the 2dHPLs, these z-dependent coefficients can contain HPLs of z, which are non-analytic at 
z — \, thus resulting in a failure of the power series expansion. In these cases, for 1/2 < z < 1, the exact 
expressions are used. 

All the various expansions (in y, but occasionally also in z) can be accelerated by the Bernoulli-like 
changes of variables jl^ already used in , and the resulting series can be further economized by standard 
use of Chebyshev polynomials. The expansion in y is always performed first, yielding z-dependent coef- 
ficients. The FORTRAN subroutine implementing the numerical evaluation of the z-dependent coefficients 
checks whether it is called repeatedly for the same value of z. In this case, the z-dependent expansion 
coefficients are not re-evaluated, which yields a considerable acceleration of running time (about a factor 
3.5). 

In the following subsections, we describe in detail the expansions used for the 2dIIPLs, depending on 
the combination of indices appearing in the 2dHPL under consideration. All these expansions have been 
generated using FORM |l^ , whose output was then converted to FORTRAN input of the required precision with 
a dedicated C program, rewriting in particular the exact coefficients generated by FORM as double precision 
floating-point numbers. 



7.1 The indices (0, 1) 

If only the indices (0) and (1) appear in m, the 2dIIPLs G{rh]y) are (up to an overall sign, which is -1-1 
for an even number of (1) in rn and —1 for an odd number of (1) in m) equal to the HPLs II(m;j/), to be 
evaluated by means of hplog |11|. 



7.2 The indices (0, 1 - z) 

If only (0) and (1 — z) appear in m(z), one can perform the mere change of variable from ti to t[ = 1^/(1 — z) 
in the multiple integral representation (2. 7), (2. 8). The individual integrands transform as follows: 



AU g(0; U) At'^ f(0; t[) , AU g(l - z; U) ^ -At'^ f(l; t'^ 



(7.2) 



and the 2dHPLs are re-expressed as HPLs of argument y/(l — z), which can then be evaluated by means 
of hplog 0. 

As an example, consider 



G(0,l-z,l-z;2/) 



diig(0;ti 

{-If 



At2 g(l - z; ia) / Ah g(l - z; ^3 



t' t' 

di'if(0;i'i) / ' At'^iil-A) r At'^i{l-t'^) 



H 0,1,1; 



y 

1 - z 



(7.3) 



7.3 The indices (0, -z) 

If only (0) and (— z) appear in to(z), the 2dHPLs can again be re-expressed as HPLs, in this case of 
argument —y/z by a mere change of variable from ti to t[ — —tij z in the multiple integral representation, 
which transforms the integrands 

dt.g(0;t.)^dt^f(0;i^), dt.g(-z;t,)^-di^f(l;i-) • (7-4) 

As a result, one obtains HPLs containing only the indices (0,1) with argument ^ — —y/z, which can then 
be evaluated by means of hplog [ pH . It is apparent that this transformation can be applied unambiguously 
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only if the trailing (or right-most) index is different from (0), which is however always the case for the 
2dHPLs in the irreducible set. 



As an example, we evaluate 

G(0, 0, -z; y) = dh g(0; h) T dt2 g(0; h) f ' dh gi-z; h) 

JQ Jo JQ 



(-1) / dt[mt[) f\t'J{Q-Q f\t'^i{l;t'^) 



JO Jo 



= -H 



(0,0,1;-^) . (7.5) 



7.4 The indices (0, 1 - z,l) 

If only the indices (0, 1 — z, 1) or (1 — z, 1) are present in m(z), then G(m(z); y) contains only the two right 
cuts in the complex plane, but no left cut. The separation of cuts is therefore not carried out. These 2dHPLs 
are expanded in s = — log(l — y/{l — z)), which corresponds to the Bernoulli-like variable appropriate to 
the cut at ?/ = (1 — z), always closer to the origin y = than the cut at ?/ = 1. As < y < (1 — z)/2, 
so that s < ln2 w 0.6931 . . ., we rescale s by a factor 4/3; in so doing, one still has (4s/3)" < 1 on the 
whole interval < y < (1 — z)/2 . The coefficients of the terms in (4s/3)", which are smaller by a factor of 
(3/4)" after rescaling, are finite-order polynomials in z, well-behaved in the whole interval < z < 1 . The 
convergence of the z-expansion is further improved by rewriting the resulting polynomials in z in terms of 
Chebyshev polynomials in {2z — 1), which finally allows coefficients smaller than the desired accuracy of 
the numerical evaluation to be identified and dropped. 
For example, the expansion of G(0, 0, 1 — z, 1; y) reads 

G(0,0,l-z,l;y) = 

3 J V 256 256 ^ / V3/ V 64 1024 ^ ^ 1024 ^ ' 

AsV / 441 1485 , , 63 , , 27 , 

Ti{2z-1)- t:^U2z - 1) - T^.n{2z 1) 



3 y V 131072 524288 ' ' 131072 ^ ' 524288 
4s \V 25749 6831 „ 3303 „ 729 



Ti(2z - 1) + — T2(2z - 1) + , .^J ^ii^z - 1) 



3 y V 65536000 26214400 ' ' 32768000 ' ' 26214400 
243 „ A /4s\^/ 7083 27 „ , 459 



T4(2z - 1) + — + ^^^^^^^^^^ + ,,,,,„,„ ri(2z - 1) - — -— — T2(2z - 1) 



65536000 ^ 7 V3/ V 419430400 83886080 ^ ' 52428800 
63 891 27 \ 

-10485760^^(2^ - 1) - 419430400^^(2^ - - 83886080^^(2z - 1) j + O [s^) (7.6) 

The presence of two non-separable right-hand cuts in functions bearing this combination of indices implies 
a somewhat slow convergence of the power series expansion. To attain the aimed double precision accuracy 
of 3 X 10~^^, one must keep 18 terms in the s-expansion of G(0, 0, 1 — z, 1; y) and expand the coefficients 
up to the 15th Chebyshev polynomial. Other functions with the same combination of indices require up to 
21 terms in the expansion, and up to the 19th Chebyshev polynomial in (2z — 1). 

7.5 The indices (0,-^,1) 

If the indices (0, — z, 1) or only (— z, 1) are present in to(z), then G(m(z);?/) contains two cuts, one left cut 
in y = — z and one right cut iny = 1. With this combination of indices, the irreducible 2dHPLs are chosen 
to always contain (1) as rightmost index. The primary cut of these functions is therefore always the right 
cut in y = 1, with y — ~z appearing only as secondary cut. 



After separating the left and right cuts according to ( |6.1| ), only G+(a; y) needs to be evaluated, since the 
secondary G_ (a; y) is always expressed in terms of known functions of lower weight; G+{d; y) is expanded in 
y dX y = and y is then replaced by the Bernoulli-type variable r = — log(l — y). The expansion coefficient 



11 



of r" is a function of z, which contains an overall factor z~" multiplying a numerator containing a finite 
order polynomial in z plus further n-th order polynomials in z times HPLs of argument z and indices of 
the subset (0, —1) only. In z 0, the coefficient and all its derivatives in z are finite. The coefficients are 
expanded in v — H{—1; z), which is the Bernoulli-type variable appropriate here. The convergence of the 
resulting series in v is further improved by re-expressing them in terms of Chebyshev polynomials, with v 
rescaled appropriately as (8w/3 — 1). 

For example, the expansion of G(0, 1, —z, 1; y) in powers of r reads 

G(0,l,-z,l;y) 

= G+(0,l,-2,l;2;) + G_(0,l,-z,l;2/) 

= + (^) (^-2zH(-l, -1; z)^ + (^)' (^+|H(-1; z) + yH(-l, -1; z)^ 

+ (LY (+—- — H(-l;z) - ^H(-l;z) - f^Hf-l, -1; z)) 
\zJ \ 18 18 ^ ^ 36 ^ ' ' 18 ^ ' '7 

/r\4 ( z^ z^ z , , 5z^ , , z"^ , \ 
+ (z) (-48 - 24 + 48^(-l^ ^) + ^^(-^^ ^) + 24^^-^= j 

/r\5/ z2 7z3 29z4 z ^ , ^ , 17^%^/ . ^ Q^z^ 
+ (z) i + 100 + 300 + 1800 - 100^(-^' ^) - W^^-^' ^) - 3600^(-^= 

+H(-l;2)G(0,l,-z;y)-2H(-l,-l;z)G(0,l;j/) . (7.7) 

The subsequent expansion of the coefficients of the r-expansion in Chebyshev polynomials in the rescaled v 
does yield rational coefficients, involving large integer numbers in numerator and denominator (recall that 
the Chebyshev expansion is performed on a power series, not on a finite order polynomial). For the sake 
of brevity, we shall only display the function truncated to order and ^ which is sufficient to highlight 
the structure of the expansion: 

G(0,l,-z,l;y) 

V 3 / V 512 128 V 3 / 512 ^ V 3 / / 

/4.y / 479^ / 8. _ X ^ /8. _ X \ 

V 3 y V 32768 8192 ^ V 3 / 32768 ^ V 3 / / 
+H(-l;z)G(0,l,-z;y)-2H(-l,-l;z)G(0,l;j/) . (7.8) 

To obtain the desired accuracy of 3 x 10^^^, one must keep 17 terms in the r-expansion of G(0, 1, — z, 1; y) 
and expand the coefficients up to the 10th Chebyshev polynomial in the rescaled v. Other functions with 
the same combination of indices require up to 18 terms in r and up to the 11th Chebyshev polynomial 
in rescaled v. Rewriting the expansion in r in terms of Chebyshev polynomials yields no improvement as 
far as the number of arithmetic operations required in the evaluation of the functions is concerned: the 
number of operations to be performed does in fact get larger in most cases. 



7.6 The indices (0, -z, 1 - z) 

If m(z) contains the indices (0, — z, 1 — z) or only (— z, 1 — z), then G(m(z);2/) contains two cuts, one left 
cut in y = — z and one right cut vtiy — 1 — z. The irreducible 2dHPLs with this combination of indices are 
chosen to contain always a (1 — z) as rightmost index. Therefore, these functions have always the right cut 
in 2/ = 1 — z as primary cut, with y — —z only appearing as secondary cut. 

After separating the left and right cuts according to (6.1), only G+(a; y) needs to be evaluated, since the 
secondary G_(a; y) is always expressed in terms of known functions of lower weight; G+(a; y) is expanded 
in 1/ at 2/ = 0, and y is then replaced by the Bernoulli-type variable s = — log(l — y/(l — z)). The 
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coefficients of the s-expansion take a form very similar to the expansion coefficients of the 2dHPLs with 
indices (0, 1 — z, 1) described in the previous subsection. The coefficient of s" is indeed equal to an overall 
factor multiplying a numerator containing a finite order polynomial plus further n-th order polynomials 
multiplying HPLs of argument z, which contain indices of the subset (0, 1) only. In z — > 0, the coefficient 
and all its derivatives in z are finite. 

In contrast to the expansion coefficients described in the previous subsection, which are analytic in the 
whole interval < z < 1, the coefficients in the s-expansion carried out here are finite, but non-analytic 
in the point z = 1. It is therefore not possible to obtain for these coefficients a series representation valid 
over the whole interval < z < 1. Instead, we cut this interval at z = 1/2. For < z < 1/2, we express 
z in terms of u = 11(1; z), which is the appropriate Bernoulli variable here, and expand the coefficients in 
u. Again, the convergence of the resulting series of powers of u is further improved by re-expressing them 
in terms of Chebyshev polynomials, with u rescaled appropriately. In the second interval 1/2 < z < 1, we 
evaluate the coefficients using the exact expressions. To avoid large-scale numerical cancellations, we then 
compute the coefficients as a function of 1 — z instead of z. 

As an example, we quote the s-expansion of G(— z. 0, 1 — z, 1 — z; y): 

G(-z,0, 1 - z,l - z;y) 

= G+(— z, 0,1 — z,l — z;y) + G_ (— z, 0, 1 — z, 1 — z; y) 

= + (^) (+zH(0, 1, 1; z) - H(0, 1, 1; z) + zH(l, 1, 1; z) - H(l, 1, 1; z)) 

+ (^) ' (+^H(0, 1, 1; z) - |h(0, 1, 1; z) + ^R{1, 1, 1; z) - ^H(l, 1, 1; z)^ 

+ ( ~ ^' ^' ^ ^' ^' " ^ ^' ^' " ^' ^' 

-f ^H(l, 1, 1; z) - ^H(l, 1, 1; z) - ^ + ^) + 0(s*) 

-HH(0, 1, 1; z)G(-z; y) + H(l, 1, 1; z)G(-z; y) . (7.9) 

Again, the subsequent expansion of the coefficients in Chebyshev polynomials in the rescaled u does yield 
rational coefficients involving large integer numbers in numerator and denominator. The structure of the 
result is very similar to that of the example quoted in the previous subsection. 

To obtain the desired accuracy of 3 x lO""*^^, one must keep 17 terms in the s expansion of G(— z, 0, 1 — 
z, 1 — z; y) and expand the coefficients up to the 10th Chebyshev polynomial in the rescaled u. Other 
functions with the same combination of indices require up to 18 terms in s and up to the 11th Chebyshev 
polynomial in the rescaled u. 



7.7 The indices (0,1, 1-2, 

The 2dHPLs develop their full analytic structure if all indices from the set (0, — z, 1 — z, 1) or (— z, 1 — z, 1) 
are present in rn(z): G(m(z); y) then contains two right cuts in y = 1 — z and y = 1, as well as one left cut 
in y = — z. The irreducible 2dHPLs with this combination of indices are chosen to contain always a (1) as 
rightmost index. Therefore, these functions have always the right cut as primary cut, with the left cut as 
secondary cut. 



Again, the left and right cuts are separated according to (6.1), and only G+{a; y) needs to be evaluated, 
since the secondary G-{d;y) is always expressed in terms of already evaluated functions of lower weight. 
G+{a; y) is expanded in y at y = 0, and y is then replaced by the Bernoulli-type variable s = — log(l — y/(l — 
z)). The expansion coefficients take a form similar to the expansion coefficients of the 2dHPLs discussed 
in the two previous subsections. The coefficient of s", which contains terms in z~", can be written as an 
overall factor z~" multiplying a numerator consisting of a finite-order polynomial in z plus further n-th 
order polynomials in z multiplying HPLs of argument z, which contain either indices of the subset (0, 1) 
only or indices of the subset (0, —1) only . In z — > 0, the coefficient and all its derivatives in z are again 
ffiiite. The simultaneous appearance of HPLs of argument z with indices (0, 1) and of HPLs with indices 
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(0, —1) forbids the introduction of either u = H{1] z) or v = H(—l; z) as Bcrnoulh variable for replacing z 
in the evaluation of the coefficients of the s-expansion. 

To evaluate those coefficients, say c„(z), in an efficient way, we write them as c„(z) = c,i.„(z) + c„_„(z), 
where c„^u{z) corresponds to the w-type contributions from HPLs of argument z and indices of the subset 
(0, 1) only, and c,j to the u-type contributions from HPLs with indices of the subset (0, —1) only. The 
separation of the contributions is evident for the polynomials in z multiplying HPLs of either combination 
of indices. The remaining polynomial in z, which does not multiply any HPL, must also be split into 
u-type and t^-type contributions, which is less evident. The underlying procedure is as follows. We recall 
that the coefficient of s" is analytic in z = 0, but has been written with an overall factor of z"". The 
remaining polynomial, as well as the terms containing HPLs of z, sum up to yield a function proportional 
to z" for z ^ 0. Expanding in z, around z = 0, the HPLs with indices (0, 1), we can identify the terms of 
the polynomial that have to attributed to c„,„(z) to ensure that c„,„(z) is analytic in z = 0. As a result, 
Cn,u{z) is finite for < z < 1; it is however non-analytic in z = 1, as the HPLs with indices (0, 1) have a 
cut at z = 1. All remaining terms of the polynomial are attributed to Cn,v{z), which is then analytic in 
the whole interval < z < 1. After performing this separation, c„,^(z) is expanded in terms of Chebyshev 
polynomials in 2z — 1, which yield an accurate description for the full interval in z. To evaluate c„^„(z), 
(which is not analytic in z = 1), we split the z-interval at z = 1/2 as in the previous subsection. Below 
this value, c„,„(z) is evaluated by replacing z by the Bernoulli- like variable u = H(l;z), and using the 
expansion in u with subsequent improvement in terms of Chebyshev polynomials. Above z = 1/2, the 
exact expression is used for the calculation. 

To illustrate the separation of u-type and u-type terms, we consider 

G(0,-z,l-z,l;y) 

= + (^) (^-^H(0, -1; z) + H(0, -1; z) + zH(0, 1; z) - H(0, 1; z)j 

+ (^) ' (+^H(0, -1; z) - hi{0, -1; z) - ^H(0, 1; z) + hi{0, 1; z)^ 

z) ( + 18 - Y + 18 + 9 12^(°' ^) - 36«(°' ^) 
-^H(0, 1; z) + ^H(0, 1; z) + |^H(0, 1; z)) + 0{s^) 
-H(0, -1; z)G(0, -z; y) + H(0, 1; z)G(0, -z; y) 
= +s (^-H(0, -1; z) + ^H(0, -1; z) j + s (^+H(0, 1; z) - ^H(0, 1; z)j 

W (+\m -1; z) - ^H(0, -1; z) + ^^+ 1; z) + ^H(0, 1; z) - 1 

W ( + 1h(0, 1; z) - ^H(0, 1; z) + J^H(0, 1; z) + ^ - ^) + 

-H(0, -1; z)G(0. -z- y) + H(0, 1; z)G(0, -z; y) 
= +SCi^y{z) + SCi^u{z) + C2^v{z) + C2,u{z) + s^C3,„(z) + s^C3,„(z) -|- 0{s'^) 

-H(0, -1; z)G(0, -z; y) -h H(0, 1; z)G(0, -z; y) . (7.10) 

An accuracy of 3 x 10""'^^ is obtained for G(0, — z, 1 — z, 1; y) if 19 terms are kept in the ,s-expansion, the 
Chebyshev expansion of the u-type contribution to the coefficients contains 18 terms, and the expansion of 
the u-type contribution another 10 terms. Other functions with the same combination of indices require up 
to 22 terms in s and up to 19 terms for the u-type contribution and 11 terms for the u-type contribution. 
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8 Checks 



We have carried out several checks of the implementation of the algorithm described in the previous section 
into the FORTRAN subroutine tdhpl. 

An immediate check of the numerical implementation of the 2dHPLs is provided by the derivative 



formula Eq. (2_^). Evaluating the left-hand side of Eq. (2^) numerically, with a standard symmetric 
4-point differentiation formula, and comparing it with the right-hand side evaluated directly, we found 
agreement within an accuracy of 10^^^ or better. This accuracy is mainly limited by rounding errors 
induced by the small interval size used in the differentiation formula. We used 10~^ as interval size, which 
implies a theoretical accuracy of about 10~^^. This accuracy is, however, reduced to the observed 10~^^ 
by rounding errors arising from taking the difference between function values evaluated at interval points 
spaced by 10~^ and evaluated only with the double precision FORTRAN accuracy. 

We also checked the continuity of all 2dHPLs across the boundaries of the two different regions intro- 
duced in Section |^, matching onto each other at y ~ (1 — z)/2. Evaluating the 2dHPLs according to the 
algorithms appropriate to the regions left and right of the boundary, we found both limiting values to agree 
within 3 x 10~^^ or better. Moreover, we evaluated the 2dHPLs in a few points scattered in a small neigh- 
bourhood (of size ±10"^°) across this boundary, by using for each point the two different algorithms (used 
separately in the FORTRAN code at each side of the boundaries) and then comparing the results. Again, we 
found agreement within the desired accuracy of 3 x 10~^^. 

Up to weight 3, all 2dHPLs can be expressed as Nielsen's generalized polylogarithms of suitable argument 
(see Appendix^). We have checked that the results produced in the above way agree numerically with their 
expressions in terms of generalized polylogarithms, evaluated using hplog |]ll|] in a number of randomly 
chosen points in y and z. 

Finally, some of the two-loop Feynman integrals computed analytically in terms of 2dHPLs in had 
been computed numerically at special points of the phase space in |l^ . We find full agreement with these 
results, within the numerical uncertainty quoted in [p^ , which is however only 1%. This comparison should 
rather be viewed as a verification of the analytical results of |0] for the two- loop Feynman integrals. 

9 The subroutine tdhpl 
9.1 Syntax 

The routine tdhpl has the following syntax: 

subroutine tdhpl (y , z , nmax , GYZl , GYZ2 , GYZ3 , GYZ4 , 
$ HZ1,HZ2,HZ3,HZ4) 



9.2 Usage 

In calling hplog, the user has to supply 

y,z: The arguments y,z for which the 2dHPLs are to be evaluated. y,z are of type real*8. They can 
take any value inside the triangle 0<?/<l — z,0<z<l- 

nmax: The maximum weight of the 2dHPLs to be evaluated, nw is of type integer. It is limited to 
1 < nw < 4. 

The output of tdhpl is provided in the arrays GYZ1,GYZ2,GYZ3,GYZ4,HZ1,HZ2,HZ3,HZ4. These have 
to be declared and dimensioned by the user as follows: 

real*8 HZ 1 , HZ2 , HZ3 , HZ4 , GYZ 1 , GYZ2 , GYZ3 , GYZ4 

dimension HZ1(0:1) ,HZ2(0 : 1 ,0 : 1) ,HZ3(0 : 1 ,0 : 1 ,0 : 1) , 
$ HZ4(0:1,0:1,0:1,0:1) 

dimension GYZl (0:3) ,GYZ2(0:3,0:3) ,GYZ3(0:3,0:3,0:3) , 
$ GYZ4(0:3,0:3,0:3,0:3) 
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It should be noted that this declaration is always needed, even if tdhpl is called with nmax< 4. After calling 
tdhpl for given arguments y,z, the arrays GYZl, GYZ2, . . . contain the values of the corresponding 
2dHPLs of weight 1,2,... , the arrays HZl , HZ2 , ... contain the values of the HPLs of argument z and 
indices (0,1), which always appear together with the 2dHPLs of argument y and index vector m(z) in 
calculations of Feynman integrals. The indices (0, 1, 1 — z, — z), which can appear in the index vector of the 
2dHPLs, correspond to the indices (0, 1,2,3) of the arrays GYZl , GYZ2 , .... 
The subroutine does not need initialization. 



9.3 Example 

The following example program illustrates how to evaluate 2dHPLs up to weight 4 for given values of y 
and z, and to write them out: 

program test2dhpl 
implicit none 
integer nmax 
integer i , il , 12 , 13 , 14 
real*8 y,z 

real*8 HZ 1 , HZ2 , HZ3 , HZ4 , GYZ 1 , GYZ2 , GYZ3 , GYZ4 
dimension HZl (0:1) ,HZ2(0 : 1 ,0 : 1) ,HZ3(0 : 1 ,0 : 1 ,0 : 1) , 
$ HZ4(0:1,0:1,0:1,0:1) , 

$ GYZl (0:3) ,GYZ2(0:3,0:3) ,GYZ3(0:3,0:3,0:3) , 

$ GYZ4(0:3,0:3,0:3,0:3) 



nmax = 4 



write (6,*) 'Input y,z:' 
read(5,*) y,z 

call tdhpl (y , z , nmax , GYZl , GYZ2 , GYZ3 , GYZ4 , HZl , HZ2 , HZ3 , HZ4) 
do il = 0,3 

write(6,101) il,GYZl(il) 
do 12 = 0,3 

write(6,102) il,i2,GYZ2(il,i2) 
do 13 = 0,3 

write (6 , 103) 11,12,13, GYZ3 (11 , 12 , 13) 
do 14 = 0,3 

write (6 , 104) il , 12 , 13 , 14 , GYZ4(il , 12 , 13 , 14) 
enddo 
enddo 
enddo 
enddo 
stop 

101 formate G(',i2,',y) = ',fl8.15) 

102 formate G(' ,12, ' , ' ,12, ' ,y) = ',fl8.15) 

103 formate GC ,12, ' , ' ,12, ' , ' ,12, ' ,y) = ',fl8.15) 

104 formate GC ,12, ' , ' ,12, ' , ' ,12, ' , ' ,12, ' ,y) = ',fl8.15) 
end 



10 Numerical examples 

In Fig. ||, we depict G(l — z, 1; y), G(0, — z, 1; y) and G(0, — z, 1 — z, 1; y) as examples of the dependence of 
the 2dHPLs on y and z. 
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G(1-z,1;y) 




Figure 2: Examples for the dependence of 2dHPLs on y and z 



11 Summary 

In this paper, we have described the routine tdhpl, which evaluates the two-dimensional harmonic poly- 
logarithms up to weight 4 for real arguments in the triangle < y < 1 — 2;, Q < z < 1. The evaluation is 
based on a series expansion in terms of appropriately transformed expansion parameters for small values 
of y < (1 — z)/2 . The evaluation for y in the interval (1 — 2:)/2<y<l— 2;is based on a transformation 
formula, relating 2dHPLs of arguments \ — y — z and y. The coefficients of the y expansion for small values 
of y, which depend on z, are then evaluated either in terms of a further power series expansion in z or using 
their exact expression. The convergence of the expansions is improved by using Bernoulli-like variables 
and Chebyshev polynomials. The algorithms used and described here can be extended to higher weights 
without further modification, requiring only the harmonic polylogarithms up to the desired weight to be 
known. 
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A Comparison of different notations 

In this appendix, we compare the different notations used for 2dHPLs in the literature. In particular, they 
form a special case of the hyperlogarithms of ^ and the multiple polylogarithms of Since we used a 
different notation for the 2dHPLs in an earlier publication [Q , we also provide the formulae to convert the 
notation used in to the notation employed here. 

A few examples for 2dIIPLs written out in the different notations are collected in Table ||. 

A.l HPLs and previously employed notation for 2dHPLs 

The 2dHPLs were introduced in |0| as a generalization of the HPLs. The one-dimensional HPL H(mu,;a;) 
is described by a lu-dimensional vector of parameters and by its argument x; w is called the weight of 
H. The HPLs are defined recursively: 



1. Definition of the HPLs ai w — I: 



H(l;a;) = -ln(l-a;), 
H(0;a;) = Inx , 
H(-l;a;) = ln(l + x); (A.l) 



definition of the rational fractions in x 

f{l;x) 
f(0;a;) 



l~x 
1 



such that 



2. For w > 1: 



X 

f(-l;x) ^ (A.2) 

1 + X 

d 

— H(a; a;) = f(a; x) with a = +1,0.-1. (A.3) 
ox 

H(0,...,0;a;) = —An^x, (A.4) 
w! 

H(a,5;a;) = [ da;'f(a; a;')H(6; x') , (A.5) 
Jo 

so that the differentiation formula is, in any case, 

— H(a, b; x) = f(a; x)Il{b; x) . (A.6) 

The notation for 2dHPLs used in 0| closely resembled the above notation of the HPLs by extending 
the set of fractions by 

f(l-z;y) : ^ 



1-y - z 

f(^;.) ^ (A.7) 
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2dHPL 


2dHPL of § 


Hyperlogarithm 


Multiple polylogarithm 


G(l - z, -z;y) 
G(0,0,-z;2/) 
G(0, -z,l-z,-z;y) 
G(l,l,l-z,-z;j/) 


-H(l - z,z;y) 
H(0,0,z;y) 
-Rio, z, I ~ z,z;y) 
-z,z;y) 


Lo{-z,l - z\y) 
Lo{-z, 1 - z, l,l\y) 


Ii,i(-z : 1 - z : y) 

hi-z ■■ y) 

li,i,2{-z : I - z : -z : y) 
: 1 - z : 1 : 1 : y) 



Table 2: Examples of 2dHPLs written in different notations. Some 2dHPLs cannot be expressed as hyper- 
logarithms. 



and correspondingly the set of HPLs at weight 1 by 

H(l-z;y) = -ln(l 



1 - 

H(z;y) = Inl^V (A.: 



Allowing (z, 1 — z) as components of the vector of parameters, ( [A.5| ) did then define the 2dHPLs. 

However, it turns out that this notation is unpleasant, since the rational fractions defined this way are 
not continuous in z as z — > 1: indeed one has 

limf(l-z;y) = -f(0;y) . (A.9) 

This nuisance is eliminated by the new notation for the 2dHPLs introduced in this paper. The relations 
between the two sets of rational factors are 

f(l;x) = -g(l;x), 

f(l - z]x) = -g(l - z;x) , 

i{z;x) = g{-z;x) , 

f(0;x) = g(0;z), (A.IO) 

and between the 2dHPLs at w = 1: 

H(l;a;) = -G(l;a;), 

H(l — z;a;) = — G(l — z;a;), 

H(z; x) = G(— z; x) , 

H(0;a;) = G(0; z) . (A.ll) 

To convert from the notation of to the new notation used here, all (z) in the index vector are to be 
replaced by (— z), and a (—1) is to be multiplied for each occurrence of (1) or (1 — z) in the index vector. 

A. 2 Hyper logarithms 

The hyperlogarithms of Lappo-Danilevsky at weight w = 1 are defined by 

\x) ^ r = log ^—^ . (A.12) 

At higher weights, the hyperlogarithms are defined recursively 

r ( I ^_ r -^fa(«Ji'°i2)---.a3.-i|a^) ^ 

Li,[aj-^,aj2, . . . ,aj^\x) — / ax. (A. Id) 
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Here & is a finite real number, which is different from any of the aj. . 

Hyper logarithms therefore include the subset of the 2dHPLs with all indices different from (0): 

G{mi{z),...,m^{z);y)^ Lo{m^{z),...,mi{z)\y) with m.,(z) 7^ . (A. 14) 

A. 3 Multiple polylogarithms 

Allowing the lower limit of integration to coincide with one or more of the elements of the index vector, 
hyperlogarithms are generalized to multiple polylogarithms. Without loss of generality, the lower limit of 
integration can be taken to be 0. The definition of multiple polylogarithms according to Goncharov Q is 
then 

Ini,. .. ,71^71 (^1 ■ ^m+l) — 

°™+^ dt dt dt dt dt dt dt dt dt 

O O . . . O O O 0...0 0...0 o — o . . . o — , (A. 15) 

t — ai t t t — a2 t t t ~ a„i t t 







mtimes n2times n,„times 

where 



dt dt f dti dt„ 







... _ _ ... . (A.16) 

The subset of the 2dHPLs with trailing index not equal to (0) is contained in the multiple polylogarithms: 

G(0„i,mi(z),0„2,m2(z), . . . , 0„^, mr(2;); y) = I(„,+i).(„,_i+i) („i+i) (TOr(^) : mr-i{z) : . . . : mi{z) : y) 

. n . . . . ^^-^S 

A similar notation has been proposed for the HPLs in |7| , and it is used in the calculations presented in [g| . 

B Relation to Nielsen's generalized polylogarithms 

For special values of the indices, it is possible to express the 2dHPLs in terms of the commonly known 
Nielsen's generalized polylogarithms Q : 



(n - 1)W in 

A special case of Sn.p{x) is the polylogarithm 

Li„(x) = S„_i,i(x) . (B.2) 

Numerical implementations of the S„_p(x) exist in the subroutine GPLDG and are widely used. 

If the index vector m(z) of G{rn{z);y) contains, besides (0), only one other index, which can appear 
more than once, but only to the right of the rightmost (0), then G{rh{z)-,y) can be expressed as 

G(0„,rp;z/) = (-lySnAv) 

y 



G(0„,^p;2/) = (-l)PS„,p (-^) (B.3) 

If, besides (0), more than one other index appears in rn[z)^ G{rn{z)]y) can be related to generalized 
polylogarithms only in special cases. Relations for all G{fh{z);y) exist only up to weight w = 3 0]. 



The irreducible 2dHPLs at w = 2 not expressed by (B.3) are as follows 



G(-M;.) ^ l„,l + *(i±i)+Lb(^)-Li,(i±i 

G{z,l-z;y) = -ln(l-z)lnf^^) +Li2(z)-Li2(?y + z) . (B.4) 
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Expressions for all irreducible 2dHPLs at w = 3 in terms of Nielsen's generalized polylogarithms are listed 
in the appendix of Q . 
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